Sesión 1
Profesor:
Dr. José
Manuel MAGALLANES REYES, Ph.D.
Profesor del Departamento de Ciencias Sociales, Sección de Ciencia Política y Gobierno.
Oficina 105 - Edificio CISEPA / ECONOMIA / CCSS
Telefono: (51) 1 - 6262000 anexo 4302
Correo Electrónico: jmagallanes@pucp.edu.pe
La Ruta hacia la Regresión
El Pipeline del Analista Político
De manera abstracta podemos ver que Los analistas existen para brindar explicaciones y recomendaciones sobre algún asunto de interés. El decisor elije el curso de acción sabiendo que no puede esperar del analista información completa. El gestor se encarga de implementar las decisiones; esa implementación traerá nueva información y el analista vuelve a su trabajo.
La estadística es una ciencia matemática que guía científicamente el análisis de datos. Esto lo hace siguiendo una secuencia:
Apuesta por entender una variable que explicaría un problema de interés (narcotráfico, elecciones, etc). En esta etapa, hace exploración de los datos de esa variable, es decir, organiza los datos en medidas, tablas y gráficos. Se entiende que la variable de interés tiene una variabilidad tal que despierta en el analista la necesidad de preguntar por qué esa variabilidad.
Se plantea hipótesis respecto a qué se relaciona con la variabilidad de la variable de interés. Las hipótesis se formula no antes de haber revisado la literatura. En esta etapa hace uso de análisis bivariado o multivariado. Esta etapa se enriquece si está actualizado en las teorías que proponen cierta asociación de la variable de interés con otras variables.
Aplica la técnica estadística que corresponda, tal que verifique la hipótesis que se planteó.
Interpretado los resultados obtenidos.
Elabora síntesis de lo actuado; propone explicaciones a lo encontrado; y elabora recomendaciones.
Hay muchas opciones para las técnicas señaladas en el punto 4. La elección de las mismas dependerá de la preparación del analista. Esta sesión te guiará para que:
I. Explorando la Variable de interés
Supongamos que estás interesado en la “situación de los locales escolares en el Perú”. Ese simple interés te lleva a buscar datos para saber tal situación. De pronto te encuentras con estos datos:
I. De Correlación a regresión
Trabajemos con estos datos:
library(rio)
organization='https://github.com/Estadistica-AnalisisPolitico'
repo='/DataFiles-estadistica/raw/main/'
file='hsb_ok.xlsx'
hsb=import(paste0(organization,repo,file))Antes de correr cualquier análisis estadístico, debes revisar como el tipo de datos que tu archivo trae es reconocido por el R:
str(hsb)## 'data.frame': 600 obs. of 15 variables:
## $ ID : chr "1" "2" "3" "4" ...
## $ SEX : num 2 1 2 2 2 1 1 2 1 2 ...
## $ RACE : num 2 2 2 2 2 2 2 2 2 2 ...
## $ SES : num 1 1 1 2 2 2 1 1 2 1 ...
## $ SCTYP : num 1 1 1 1 1 1 1 1 1 1 ...
## $ HSP : num 3 2 2 3 3 2 1 1 1 1 ...
## $ LOCUS : num 0.29 -0.42 0.71 0.06 0.22 0.46 0.44 0.68 0.06 0.05 ...
## $ CONCPT: num 0.88 0.03 0.03 0.03 -0.28 0.03 -0.47 0.25 0.56 0.15 ...
## $ MOT : num 0.67 0.33 0.67 0 0 0 0.33 1 0.33 1 ...
## $ CAR : num 10 2 9 15 1 11 10 9 9 11 ...
## $ RDG : num 33.6 46.9 41.6 38.9 36.3 49.5 62.7 44.2 46.9 44.2 ...
## $ WRTG : num 43.7 35.9 59.3 41.1 48.9 46.3 64.5 51.5 41.1 49.5 ...
## $ MATH : num 40.2 41.9 41.9 32.7 39.5 46.2 48 36.9 45.3 40.5 ...
## $ SCI : num 39 36.3 44.4 41.7 41.7 41.7 63.4 49.8 47.1 39 ...
## $ CIV : num 40.6 45.6 45.6 40.6 45.6 35.6 55.6 55.6 55.6 50.6 ...Casi todo salió numérico, pero no sabremos qué ajustar si no leemos el manual metodológico o el diccionario de datos o la metadata disponible.
De ahi que debemos pre procesar:
categoricals=c("SEX","RACE","SES","SCTYP","HSP","CAR")
hsb[,categoricals]=lapply(hsb[,categoricals], as.factor)
# nominales
hsb$SEX=factor(hsb$SEX,
levels=c(1,2),
labels=c("Male","Female"))
hsb$RACE=factor(hsb$RACE,
levels=c(1,2,3,4),
labels=c("Hispanic","Asian","Black","White"))
hsb$HSP=factor(hsb$HSP,
levels=c(1,2,3),labels=c("General","Academic","Vocational"))
hsb$SCTYP=factor(hsb$SCTYP,
levels=c(1,2),
labels=c("Public","Private"))
# a ordinal:
hsb$SES=ordered(hsb$SES,
levels=c(1,2,3),
labels=c("Low","Medium","High" ))La variable de interés:
Asumamos que nuestra variable de interés es el desempeño en matemáticas; así, nuestra variable dependiente está representada por la variable MATH.
A partir de ahí, consideremos que nos interesa saber la posible relación que pueda tener la variable que ha medido el desempeño en escritura; así, una variable independiente sería la representada por la variable WRTG. Hasta ahora sabemos que como son dos variables de tipo numérico debemos usar una correlación. La gráfica de correlación es esta:
library(ggplot2)
base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
scatter = base + geom_point()
scatter
Vemos que hay una aparente relación. Calculemos los indices de correlación:
f1=formula(~MATH + WRTG)
# camino parametrico
pearsonf1=cor.test(f1,data=hsb)[c('estimate','p.value')]
# camino no parametrico
spearmanf1=cor.test(f1,data=hsb,method='spearman')[c('estimate','p.value')]Asumiendo un camino paramétrico, podemos pedir el coeficiente de Pearson, el cuál al ser calculado obtenemos 0.6326664 (con p-value= 0).
Si hubieramos seguido una ruta no paramétrica, informaríamos el coeficiente de Spearman:0.6415126 (con p-value= 0).
Ahora, consideremos que nos interesa saber a la vez la posible relación que pueda tener la variable que ha medido el desempeño en ciencias; así otra variable independiente sería la representada por la variable SCI. Como es otra variable numérica no podemos calcular la correlación de tres variables, pero podemos tratar de verlo visualmente:
- En tres dimensiones:
library(scatterplot3d)
scatterplot3d(hsb[,c('SCI','WRTG','MATH')])
- En dos dimensiones:
base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
base + geom_point(aes(color = SCI))
Calculemos las correlaciones:
f2=formula(~MATH+SCI)
# camino parametrico
pearsonf2=cor.test(f2,data=hsb)[c('estimate','p.value')]
# camino no parametrico
spearmanf2=cor.test(f2,data=hsb, method='spearman')[c('estimate','p.value')]Podríamos calcular la correlación de SCI con MATH, obteniendo el Pearson (0.6495261, p-value= 0) y el Spearman (0.6551515,p-value= 0).
Visualmente vemos relación, pero no tenemos un coeficiente para medir ello.
Finalmente, ¿si quisiéramos ver si el sexo tiene algun rol en todo esto? Como ésta es una variable categórica y dicotómica, lo primero que puede venir a la mente es esta gráfica:
base=ggplot(data=hsb, aes(x=as.factor(SEX), y=MATH))
base + geom_boxplot(notch = T) + geom_jitter(color="black", size=0.4, alpha=0.9)
Los boxplots tienen un notch flanqueando a la mediana, para sugerir igualdad de medianas si éstos se intersectan; de ahi que parece no haber diferencia sustantiva entre hombres y mujeres en cuanto a su desempeño en MATH.
Una alternativa al boxplot seria las barras de error:
library(ggpubr)
ggerrorplot(data=hsb, x = "SEX", y = "MATH")
En este último caso, si las lineas (denotado por las barras de error de la media) se intersectan, sugeriria que los valores medios (en este caso la media) podrian ser iguales.
Verificar si hay o no igualdad entre distribuciones depende si las variables se distribuyen o no de manera normal:
library(ggplot2)
ggplot(hsb,aes(x=MATH)) +
geom_histogram(aes(y = ..density..),bins = 20, fill='green') +
stat_function(fun = dnorm, colour = "red",
args = list(mean = mean(hsb$MATH, na.rm = TRUE),
sd = sd(hsb$MATH, na.rm = TRUE))) +
facet_grid(~SEX) +
coord_flip()
Nota que los histogramas de la data real tienen encima la curva normal que idealmente tendría esa data. La lejanía entre ellos, sugeriría no normalidad.
Se suele usar un qqplot para explorar la presencia/ausencia de normalidad:
# se sugiere normalidad si los puntos no se alejan de la diagonal.
library(ggpubr)
ggqqplot(data=hsb,x="MATH") + facet_grid(. ~ SEX)
Como ello no es fácil de discernir visualmente, tenemos por costumbre calcular algun coeficiente, como el Shapiro-Wilk:
library(knitr)
library(magrittr)
library(kableExtra)
f4=formula(MATH~SEX)
tablag= aggregate(f4, hsb,
FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})
# para que se vea mejor:
shapiroTest=as.data.frame(tablag[,2])
names(shapiroTest)=c("W","Prob")
kable(cbind(tablag[1],shapiroTest))%>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = F, position = "left")| SEX | W | Prob |
|---|---|---|
| Male | 0.9837903 | 0.0034565 |
| Female | 0.9790040 | 0.0001031 |
Esto nos sugiere un camino no paramétrico para ver la diferencia de valores medios, por lo que deberiamos usar la prueba de Mann-Whitney en vez de la prueba t para testaer la relación entre ambas.
tf4=t.test(f4,data=hsb)[c('estimate','p.value')]
wilcoxf4=wilcox.test(f4,data=hsb)['p.value']La prueba no paramétrica no rechazaría la igualdad de valores medios (Mann-Whitney con p valor = 0.3085543).
- Veamos como representar al sexo en nuestra gráfica entre WRTG y MATH:
base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
base + geom_point(aes(color = SEX))
- Veamos como representar al sexo en nuestra gráfica entre WRTG, SCI y MATH:
base + geom_point(aes(size = SCI, color=SEX)) 
Otra alternativa puede ser:
base + geom_point(aes(color = SCI)) + facet_grid(~SEX)
Y claro:
paleta <- c("coral1","cyan" )
colors <- paleta[as.numeric(hsb$SEX)]
scatterplot3d(hsb[,c('SCI','WRTG','MATH')],color=colors)
En todos los gráficos vemos que los hombres y las mujeres están distribuidos por todo el gráfico, lo cual nos sugiere que no hay diferencias aun en dimensiones mayores a dos. Sin embargo, no tenemos una medida de cuanto cada uno afecta a nuestra dependiente.
De ahi que necesitamos la regresión.
Regresión Lineal
La regresión es una técnica donde hay que definir una variable dependiente y una o más independientes. Las independientes pueden tener rol predictor, dependiendo del diseño de investigación, aunque por defecto tiene un rol explicativo.
La regresión sí quiere informar cuánto una variable (independiente) puede explicar la variación de otra (dependiente), de ahí que es una técnica para probar hipótesis direccionales o asimétricas (las correlaciones tiene hipótesis simétricas).
La regresión busca proponer un modelo, es decir una ecuación, que recoja como una variable explicaría a otra.Para nuestro caso la variable dependiente es MATH, y tendriamos hasta ahora tres modelos:
modelo1=formula(MATH~WRTG)
modelo2=formula(MATH ~ WRTG + SCI)
modelo3= formula(MATH ~ WRTG + SCI + SEX)Por ejemplo, para la hipótesis ‘el nivel de desempeño en escritura afecta el desempeño en matemáticas’, la regresión arrojaría este resultado:
O en mejor version con ayuda de stargazer:
library(stargazer)
reg1=lm(modelo1,data=hsb)
stargazer(reg1,type = "html",intercept.bottom = FALSE)| Dependent variable: | |
| MATH | |
| Constant | 19.769*** |
| (1.633) | |
| WRTG | 0.612*** |
| (0.031) | |
| Observations | 600 |
| R2 | 0.400 |
| Adjusted R2 | 0.399 |
| Residual Std. Error | 7.297 (df = 598) |
| F Statistic | 399.110*** (df = 1; 598) |
| Note: | p<0.1; p<0.05; p<0.01 |
Aquí ya sabemos algo interesante, primero que WRTG tiene efecto, pues es significativo (indicado por los dos asteriscos); segundo, que ese efecto es directo, pues el coeficiente calculado es positivo; y tercero que la magnitud de ese efecto es 0.612, lo que indica cuanto aumenta MATH en promedio cuando WRTG se incremente en una unidad.
Esto es información suficiente para representar esa relación con una ecuación. Como la ecuación sólo tiene una variable independiente podemos producir una recta sobre el gráfico de correlación:
ggplot(hsb, aes(x=WRTG, y=MATH)) +
geom_point()+
geom_smooth(method=lm)
Esa recta podemos representarla así:
\[ MATH= 19.7690323 + 0.6123904 \cdot WRTG + \epsilon\]
El Y verdadero es MATH, pero la regresión produce un \(\hat{MATH}\) estimado, de ahi la presencia del \(\epsilon\). Justamente el R cuadrado ajustado (0.4002668) nos brinda un porcentaje (multiplicalo por 100) que da una pista de nuestra cercanía a una situación perfecta (cuando vale 1).
Y sí queremos ver el efecto de SCI?
reg2=lm(modelo2,data=hsb)
stargazer(reg2,type = "html",intercept.bottom = FALSE)| Dependent variable: | |
| MATH | |
| Constant | 10.629*** |
| (1.630) | |
| WRTG | 0.377*** |
| (0.033) | |
| SCI | 0.415*** |
| (0.033) | |
| Observations | 600 |
| R2 | 0.524 |
| Adjusted R2 | 0.523 |
| Residual Std. Error | 6.505 (df = 597) |
| F Statistic | 328.846*** (df = 2; 597) |
| Note: | p<0.1; p<0.05; p<0.01 |
En este caso, la regresión tendrá una formula con dos variables explicando la dependiente, así que en vez de producir una línea buscará producir un plano:
library(scatterplot3d)
G <- scatterplot3d(hsb[,c('SCI','WRTG','MATH')])
G$plane3d(reg2, draw_polygon = TRUE, draw_lines = FALSE)
Este plano podemos representarlo así:
\[ MATH= 10.6285904 + 0.3765304 \cdot WRTG + 0.4152733 \cdot SCI + \epsilon\]
Nuevamente, el Y verdadero es MATH, pero la regresión produce un \(\hat{MATH}\) estimado en forma de plano. De igual manera el R cuadrado ajustado (0.5241861) nos da una pista de nuestra lejanía a una situación perfecta.
Es clave darse cuenta de otro detalle, que el coeficiente de WRTG ha variado en la fórmula ahora que está presente SCI ¿Por qué sucede esto? Veamoslo así: en el primer caso, WRTG y \(\epsilon\) buscaban representar la variabilidad en MATH, y ahora, en el segundo caso, viene SCI para mejorar esa explicación; es así que el peso de la explicación ahora se recalcula y el coeficiente de WRTG deja de explicar lo que le corresponde a SCI, y \(\epsilon\) también le entrega algo a SCI.
Como \(\epsilon\) no tiene coeficiente, representamos su variación usando el error típico de los residuos o residual standard error (RSE). Nótese que éste ha variado de un modelo ha otro, ahora tenemos un RSE menor. Aquí vale la pena preguntarse si esta disminución del error es significativa, obteniendo:
tanova=anova(reg1,reg2)
stargazer(tanova,type = 'html',summary = F,title = "Table de Análisis de Varianza")| Res.Df | RSS | Df | Sum of Sq | F | Pr(> F) | |
| 1 | 598 | 31,842.070 | ||||
| 2 | 597 | 25,262.730 | 1 | 6,579.335 | 155.481 | 0 |
La comparación de modelos usando la tabla de análisis de varianza (anova) propone como hipótesis nula que los modelos no difieren (no se ha reducido el error al pasar de un modelo a otro). Como la comparación es significativa (vea el P), rechazamos igualdad de modelos: el modelo 2 sí reduce el error al incluir una variable más.
Finalmente, veamos el rol de sexo:
reg3=lm(modelo3,data=hsb)
stargazer(reg3,type = "html",intercept.bottom = FALSE)| Dependent variable: | |
| MATH | |
| Constant | 11.306*** |
| (1.630) | |
| WRTG | 0.424*** |
| (0.036) | |
| SCI | 0.375*** |
| (0.035) | |
| SEXFemale | -1.922*** |
| (0.582) | |
| Observations | 600 |
| R2 | 0.533 |
| Adjusted R2 | 0.530 |
| Residual Std. Error | 6.452 (df = 596) |
| F Statistic | 226.510*** (df = 3; 596) |
| Note: | p<0.1; p<0.05; p<0.01 |
Aunque no podemos graficar cuatro coordenadas, podemos usar elementos visuales:
library(scatterplot3d)
colors <- paleta[as.numeric(hsb$SEX)]
G <- scatterplot3d(hsb[,c('SCI','WRTG','MATH')],color=colors)
G$plane3d(reg2, draw_polygon = TRUE, draw_lines = FALSE)
Nuestra nueva ecuación sería:
\[ MATH= 11.306349 + 0.4235786 \cdot WRTG + 0.3748025 \cdot SCI + -1.9219692 \cdot SEX + \epsilon\]
Nuevamente podemos ver si añadir SEXO en este modelo representa una mejora al anterior:
tanova=anova(reg1,reg2,reg3)
stargazer(tanova,type = 'html',summary = F,title = "Table de Análisis de Varianza")| Res.Df | RSS | Df | Sum of Sq | F | Pr(> F) | |
| 1 | 598 | 31,842.070 | ||||
| 2 | 597 | 25,262.730 | 1 | 6,579.335 | 158.063 | 0 |
| 3 | 596 | 24,808.420 | 1 | 454.308 | 10.914 | 0.001 |
Finalmente, podemos resumir todo en esta tabla:
library(stargazer)
stargazer(reg1,reg2,reg3, type = "html", title = "Modelos planteadas",digits = 2, single.row = F,no.space = F,intercept.bottom = FALSE,
dep.var.caption="Variable dependiente:",
dep.var.labels="Desempeño en Matemáticas",
covariate.labels=c("Constante","Desempeño en Escritura","Desempeño en Ciencias","SEXO (mujer)"),
keep.stat = c("n","adj.rsq","ser"),df = F,
notes.label = "Notas:")| Variable dependiente: | |||
| Desempeño en Matemáticas | |||
| (1) | (2) | (3) | |
| Constante | 19.77*** | 10.63*** | 11.31*** |
| (1.63) | (1.63) | (1.63) | |
| Desempeño en Escritura | 0.61*** | 0.38*** | 0.42*** |
| (0.03) | (0.03) | (0.04) | |
| Desempeño en Ciencias | 0.42*** | 0.37*** | |
| (0.03) | (0.04) | ||
| SEXO (mujer) | -1.92*** | ||
| (0.58) | |||
| Observations | 600 | 600 | 600 |
| Adjusted R2 | 0.40 | 0.52 | 0.53 |
| Residual Std. Error | 7.30 | 6.51 | 6.45 |
| Notas: | p<0.1; p<0.05; p<0.01 | ||
Gráficamente:
library(ggplot2)
library(sjPlot)## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!plot_models(reg1,reg2,reg3,vline.color = "grey",m.labels=c("Modelo 1","Modelo 2","Modelo 3"))
Como vemos, el propósito del último gráfico es mostrar que dado que ninguna de las variables (y sus intervalos de confianza) tocan el valor cero (que significa que la variable no tiene efecto en la dependiente).